/* Copyright 2007 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 * Marcus Cuda - http://marcus.cuda.net
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;
using System.Collections.Generic;
using dnAnalytics.LinearAlgebra.Decomposition;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra
{
    /// <summary>
    /// A matrix class using dense storage.
    /// </summary>
    public partial class DenseMatrix : Matrix
    {
        private readonly double[] mData;

        /// <summary>
        /// Creates a matrix with the given <paramref name="order"/>. The matrix values are the set to zero.
        /// </summary>
        /// <param name="order">The order of the new matrix.</param>
        /// <exception cref="ArgumentException">If either <paramref name="order"/> is less than one.</exception>
        public DenseMatrix(int order) : this(order, order)
        {
        }

        /// <summary>
        /// Creates a  <paramref name="rows"/> by <paramref name="columns"/> matrix. The matrix values are the set to zero.
        /// </summary>
        /// <param name="rows">The number of rows of the new matrix.</param>
        /// <param name="columns">The number of columns of the new matrix.</param>
        /// <exception cref="ArgumentException">If either <paramref name="rows"/> or <paramref name="columns"/>
        /// is less than one.</exception>
        public DenseMatrix(int rows, int columns) : base(rows, columns)
        {
            mData = new double[rows*columns];
        }

        /// <summary>
        /// Creates a rows by columns <strong>Matrix</strong>. The matrix values are the set to the given
        /// value.
        /// </summary>
        /// <param name="rows">The number of rows of the new matrix.</param>
        /// <param name="columns">The number of columns of the new matrix.</param>
        /// <param name="value">The value to set all elements of the matrix to.</param>
        /// <exception cref="ArgumentException">If either <paramref name="rows"/> or <paramref name="columns"/>
        /// is less than one.</exception>
        public DenseMatrix(int rows, int columns, double value) : this(rows, columns)
        {
            for (int i = 0; i < mData.Length; i++)
            {
                mData[i] = value;
            }
        }

        /// <summary>
        /// Creates a new matrix from another matrix.
        /// </summary>
        /// <param name="other">The matrirx to create the new matrix from.</param>
        public DenseMatrix(Matrix other) : this(other.Rows, other.Columns)
        {
            DenseMatrix matrix = other as DenseMatrix;
            if (matrix == null)
            {
                //using enumerators since they will be more efficient for copying sparse matrices
                foreach (KeyValuePair<int, Vector> column in other.GetColumnEnumerator())
                {
                    int jIndex = column.Key*Rows;
                    foreach (KeyValuePair<int, double> element in column.Value)
                    {
                        mData[jIndex + element.Key] = element.Value;
                    }
                }
            }
            else
            {
                Buffer.BlockCopy(matrix.mData, 0, mData, 0, mData.Length*Constants.SizeOfDouble);
            }
        }


        /// <summary>
        /// Creates a matrix from a multidimensional array.
        /// </summary>
        /// <param name="array">The arrau to create the matrix from.</param>
        public DenseMatrix(double[,] array) : this(array.GetLength(0), array.GetLength(1))
        {
            for (int j = 0; j < Columns; j++)
            {
                int jIndex = j*Rows;
                for (int i = 0; i < Rows; i++)
                {
                    mData[jIndex + i] = array[i, j];
                }
            }
        }

        internal double[] Data
        {
            get { return mData; }
        }

        /// <summary>
        /// Retrieves the value at the given indices. <strong>This method is not ranged checked.</strong> If the
        /// indices are outside the dimensions of the matrix, it may return an incorrect value or throw
        /// an <see cref="IndexOutOfRangeException"/>.
        /// </summary>
        /// <param name="row">The row index.</param>
        /// <param name="column">The column index.</param>
        /// <returns>The value at the given indices.</returns>
        /// <exception cref="IndexOutOfRangeException">Depending on the implementation, an <see cref="IndexOutOfRangeException"/>
        /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
        protected internal override double ValueAt(int row, int column)
        {
            return mData[column*Rows + row];
        }

        /// <summary>
        /// Set the value at the given indices to the given value. <strong>This method is not ranged checked.</strong> If the
        /// indices are outside the dimensions of the matrix, it may set an incorrect element or throw
        /// an <see cref="IndexOutOfRangeException"/>.
        /// </summary>
        /// <param name="row">The row index.</param>
        /// <param name="column">The column index.</param>
        /// <param name="value">The value to set.</param>
        /// <exception cref="IndexOutOfRangeException">Depending on the implementation, an <see cref="IndexOutOfRangeException"/>
        /// may be thrown if one of the indices is outside the dimensions of the matrix.</exception>
        protected internal override void ValueAt(int row, int column, double value)
        {
            mData[column*Rows + row] = value;
        }

        /// <summary>
        /// Creates a <strong>Matrix</strong> for the given number of rows and columns.
        /// </summary>
        /// <param name="numberOfRows">The number of rows.</param>
        /// <param name="numberOfColumns">The number of columns.</param>
        /// <returns>
        /// A <strong>Matrix</strong> with the given dimensions.
        /// </returns>
        protected internal override Matrix CreateMatrix(int numberOfRows, int numberOfColumns)
        {
            return new DenseMatrix(numberOfRows, numberOfColumns);
        }

        /// <summary>
        /// Creates a <strong>Vector</strong> with a the given dimension.
        /// </summary>
        /// <param name="size">The size of the vector.</param>
        /// <returns>
        /// A <strong>Vector</strong> with the given dimension.
        /// </returns>
        protected internal override Vector CreateVector(int size)
        {
            return new DenseVector(size);
        }

        /// <summary>
        /// Returns a deep-copy clone of the matrix.
        /// </summary>
        /// <returns>A deep-copy clone of the matrix.</returns>
        public override Matrix Clone()
        {
            return new DenseMatrix(this);
        }

        /// <summary>
        /// Sets all values to zero.
        /// </summary>
        public override void Clear()
        {
            Array.Clear(mData, 0, mData.Length);
        }

        /// <summary>
        /// Copies the elements of this matrix to the given matrix.
        /// </summary>
        /// <param name="target">The matrix to copy values into.</param>
        /// <exception cref="ArgumentNullException">If target is <strong>null</strong>.</exception>
        /// <exception cref="NotConformableException">If target is not the same size as this matrix.</exception>
        public override void CopyTo(Matrix target)
        {
            if (target == null)
            {
                throw new ArgumentNullException("target");
            }
            if (Rows != target.Rows || Columns != target.Columns)
            {
                throw new NotConformableException();
            }

            DenseMatrix dense = target as DenseMatrix;
            if (dense == null)
            {
                base.CopyTo(target);
            }
            else
            {
                Buffer.BlockCopy(mData, 0, dense.mData, 0, mData.Length*Constants.SizeOfDouble);
            }
        }

        /// <summary>
        /// Returns the transpose of this matrix.
        /// </summary>        
        /// <returns>The transpose of this matrix.</returns>
        public override Matrix Transpose()
        {
            DenseMatrix ret = new DenseMatrix(Columns, Rows);
            for (int j = 0; j < Columns; j++)
            {
                int index = j*Rows;
                for (int i = 0; i < Rows; i++)
                {
                    ret.mData[i*Columns + j] = mData[index + i];
                }
            }
            return ret;
        }

        /// <summary>
        /// Returns the matrix's elements as an array with the data laid out column-wise.
        /// </summary>
        /// <example><pre>
        /// 1, 2, 3
        /// 4, 5, 6  will be returned as  1, 4, 7, 2, 5, 8, 3, 6, 9
        /// 7, 8, 9
        /// </pre></example>
        /// <returns>An array containing the matrix's elements.</returns>
        public override double[] ToColumnWiseArray()
        {
            double[] ret = new double[mData.Length];
            Buffer.BlockCopy(mData, 0, ret, 0, mData.Length*Constants.SizeOfDouble);
            return ret;
        }


        /// <summary>
        /// Copies the values of the given array to the specified column.
        /// </summary>
        /// <param name="index">The column to copy the values to.</param>
        /// <param name="source">The array to copy the values from.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="source"/> is <see langword="null" />.</exception>        
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="index"/> is less than zero,
        /// or greater than or equal to the number of columns.</exception>
        /// <exception cref="NotConformableException">If the size of <paramref name="source"/> does not
        /// equal the number of rows of this <strong>Matrix</strong>.</exception>
        /// <exception cref="NotConformableException">If the size of <paramref name="source"/> does not
        /// equal the number of rows of this <strong>Matrix</strong>.</exception>
        public override void SetColumn(int index, double[] source)
        {
            if (index < 0 || index >= Columns)
            {
                throw new ArgumentOutOfRangeException("index");
            }

            if (source == null)
            {
                throw new ArgumentNullException("source");
            }

            if (source.Length != Rows)
            {
                throw new NotConformableException("source", Resources.ParameterNotConformable);
            }

            Buffer.BlockCopy(source, 0, mData, index*Rows*Constants.SizeOfDouble, source.Length*Constants.SizeOfDouble);
        }

        /// <summary>
        /// Copies the requested column elements into the given vector.
        /// </summary>
        /// <param name="columnIndex">The column to copy elements from.</param>
        /// <param name="rowIndex">The row to start copying from.</param>
        /// <param name="length">The number of elements to copy.</param>
        /// <param name="result">The <see cref="Vector"/> to copy the column into.</param>
        /// <exception cref="ArgumentNullException">If the result <see cref="Vector"/> is <see langword="null" />.</exception>
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="columnIndex"/> is negative,
        /// or greater than or equal to the number of columns.</exception>        
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> is negative,
        /// or greater than or equal to the number of rows.</exception>        
        /// <exception cref="ArgumentOutOfRangeException">If <paramref name="rowIndex"/> + <paramref name="length"/>  
        /// is greater than or equal to the number of rows.</exception>
        /// <exception cref="ArgumentException">If <paramref name="length"/> is not positive.</exception>
        /// <exception cref="NotConformableException">If <strong>result.Count lt; length</strong>.</exception>
        public override void GetColumn(int columnIndex, int rowIndex, int length, Vector result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (columnIndex >= Columns || columnIndex < 0)
            {
                throw new ArgumentOutOfRangeException("columnIndex");
            }
            if (rowIndex >= Rows || rowIndex < 0)
            {
                throw new ArgumentOutOfRangeException("rowIndex");
            }

            if (rowIndex + length > Rows)
            {
                throw new ArgumentOutOfRangeException("length");
            }

            if (length < 1)
            {
                throw new ArgumentException(Resources.NotPositive, "length");
            }

            if (result.Count < length)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if( result is DenseVector)
            {
                Buffer.BlockCopy(Data, (columnIndex*Rows + rowIndex)*Constants.SizeOfDouble, ((DenseVector)result).Data, 0, length*Constants.SizeOfDouble);                
            }else
            {
                base.GetColumn(columnIndex, rowIndex, length, result);
            }
        }

        /// <summary>
        /// Concatenates this matrix with the given matrix and places the result into the result matrix.
        /// </summary>
        /// <param name="right">The matrix to concatenate.</param>
        /// <param name="result">The combined matrix.</param>
        /// <exception cref="ArgumentNullException">If right is <see langword="null"/>.</exception>
        /// <exception cref="ArgumentNullException">If the result matrix is <see langword="null"/>.</exception>
        /// <exception cref="NotConformableException">If <strong>this.Rows != right.Rows</strong>.</exception>
        /// <exception cref="NotConformableException">If the result matrix's dimensions are not this.Rows x (this.Columns + right.Columns).</exception>
        public override void Append(Matrix right, Matrix result)
        {
            if (right == null)
            {
                throw new ArgumentNullException("right");
            }

            if (right.Rows != Rows)
            {
                throw new NotConformableException("right", Resources.ParameterNotConformable);
            }

            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (result.Columns != (Columns + right.Columns) || result.Rows != Rows)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            int sizeThis = Rows*Columns*Constants.SizeOfDouble;
            int sizeRight = right.Rows*right.Columns*Constants.SizeOfDouble;
            DenseMatrix denseResult = (DenseMatrix) result;
            Buffer.BlockCopy(mData, 0, denseResult.mData, 0, sizeThis);
            Buffer.BlockCopy(((DenseMatrix)right).mData, 0, denseResult.mData, sizeThis, sizeRight);
        }

        /// <summary>
        /// Adds a scalar to each element in the matrix.
        /// </summary>
        /// <param name="scalar">The scalar to add.</param>        
        public override void Add(double scalar)
        {
            if (FloatingPointUtilities.EqualsWithinSignificantFigures(0.0, scalar, 15))
            {
                return;
            }

            for (int i = 0; i < mData.Length; i++)
            {
                mData[i] += scalar;
            }
        }

        /// <summary>Computes the determinant of this matrix.</summary>
        /// <returns>The determinant of this matrix.</returns>
        /// <exception cref="MatrixNotSquareException">If the matrix is not square.</exception>
        /// <remarks>The determinant is calculated using LU decomposition.</remarks>
        public override double Determinant()
        {
            if (Rows != Columns)
            {
                throw new MatrixNotSquareException();
            }
            DenseLU lu = new DenseLU(this);
            return lu.Determinant();
        }

        /// <summary>
        /// Returns the inverse of this matrix. The inverse is calculated using LU decomposition.
        /// </summary>
        /// <returns>The inverse of this matrix.</returns>
        /// <exception cref="MatrixNotSquareException">If this matrix is not square.</exception>
        public override Matrix Inverse()
        {
            if (Rows != Columns)
            {
                throw new MatrixNotSquareException();
            }
            return new DenseLU(this).Inverse();
        }

        /// <summary>
        /// Adds two matrices together and returns the results.
        /// </summary>
        /// <param name="leftSide">One of the matrices to add.</param>
        /// <param name="rightSide">One of the matrices to add.</param>
        /// <returns>The result of the addition.</returns>
        /// <exception cref="NotConformableException">If <paramref name="leftSide"/> and <paramref name="rightSide"/> don't have the same dimensions.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
        public static DenseMatrix operator +(DenseMatrix leftSide, DenseMatrix rightSide)
        {
            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }

            if (leftSide.Rows != rightSide.Rows || leftSide.Columns != rightSide.Columns)
            {
                throw new NotConformableException(Resources.ParametersNotConformable);
            }

            DenseMatrix ret = new DenseMatrix(leftSide);
            ret.Add(rightSide);
            return ret;
        }

        /// <summary>
        /// Negates each element of the matrix.
        /// </summary>
        /// <param name="rightSide">The matrix to negate.</param>
        /// <returns>A matrix containing the negated values.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="rightSide"/> is <see langword="null" />.</exception>
        public static DenseMatrix operator -(DenseMatrix rightSide)
        {
            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }
            DenseMatrix ret = new DenseMatrix(rightSide);
            ret.Negate();
            return ret;
        }

        /// <summary>
        /// Subtracts two matrices and returns the results.
        /// </summary>
        /// <param name="leftSide">The matrix to subtract from.</param>
        /// <param name="rightSide">The matrix to subtract.</param>
        /// <returns>The result of the subtraction.</returns>
        /// <exception cref="NotConformableException">If <paramref name="leftSide"/> and <paramref name="rightSide"/> don't have the same dimensions.</exception>
        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
        public static DenseMatrix operator -(DenseMatrix leftSide, DenseMatrix rightSide)
        {
            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }

            if (leftSide.Rows != rightSide.Rows || leftSide.Columns != rightSide.Columns)
            {
                throw new NotConformableException(Resources.ParametersNotConformable);
            }

            DenseMatrix ret = new DenseMatrix(leftSide);
            ret.Subtract(rightSide);
            return ret;
        }

        /// <summary>
        /// Divides a <strong>Matrix</strong> by a constant and returns the result.
        /// </summary>
        /// <param name="leftSide">The matrix to divide.</param>
        /// <param name="rightSide">The constant to divide the matrix by.</param>
        /// <returns>The result of the division.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> is <see langword="null" />.</exception>
        public static DenseMatrix operator /(DenseMatrix leftSide, double rightSide)
        {
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }
            DenseMatrix ret = new DenseMatrix(leftSide);
            ret.Divide(rightSide);
            return ret;
        }

        /// <summary>
        /// Multiplies a <strong>Matrix</strong> by a constant and returns the result.
        /// </summary>
        /// <param name="leftSide">The matrix to multiply.</param>
        /// <param name="rightSide">The constant to multiply the matrix by.</param>
        /// <returns>The result of the multiplication.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> is <see langword="null" />.</exception>
        public static DenseMatrix operator *(DenseMatrix leftSide, double rightSide)
        {
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }
            DenseMatrix ret = new DenseMatrix(leftSide);
            ret.Multiply(rightSide);
            return ret;
        }

        /// <summary>
        /// Multiplies a <strong>Matrix</strong> by a constant and returns the result.
        /// </summary>
        /// <param name="leftSide">The matrix to multiply.</param>
        /// <param name="rightSide">The constant to multiply the matrix by.</param>
        /// <returns>The result of the multiplication.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="rightSide"/> is <see langword="null" />.</exception>
        public static DenseMatrix operator *(double leftSide, DenseMatrix rightSide)
        {
            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }
            DenseMatrix ret = new DenseMatrix(rightSide);
            ret.Multiply(leftSide);
            return ret;
        }

        /// <summary>
        /// Multiplies two matrices.
        /// </summary>
        /// <param name="leftSide">One of the matrices to multiply.</param>
        /// <param name="rightSide">One of the matrices to multiply.</param>
        /// <returns>The result of multiplication.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If the dimensions of <paramref name="leftSide"/> or <paramref name="rightSide"/> don't conform.</exception>
        public static DenseMatrix operator *(DenseMatrix leftSide, DenseMatrix rightSide)
        {
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }

            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }

            if (leftSide.Columns != rightSide.Rows)
            {
                throw new NotConformableException(Resources.ParametersNotConformable);
            }
            return (DenseMatrix)leftSide.Multiply(rightSide);
        }

        /// <summary>
        /// Multiplies a <strong>Matrix</strong> and a <see cref="Vector"/>.
        /// </summary>
        /// <param name="leftSide">The matrix to multiply.</param>
        /// <param name="rightSide">The vector to multiply.</param>
        /// <returns>The result of multiplication.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If <strong>leftSide.Columns ! = rightSide.Count</strong>.</exception>
        public static DenseVector operator *(DenseMatrix leftSide, DenseVector rightSide)
        {
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }

            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }

            if (leftSide.Columns != rightSide.Count)
            {
                throw new NotConformableException(Resources.ParametersNotConformable);
            }

            return (DenseVector)leftSide.Multiply(rightSide);
        }

        /// <summary>
        /// Multiplies a <see cref="Vector"/> and a <strong>Matrix</strong>.
        /// </summary>
        /// <param name="leftSide">The vector to multiply.</param>
        /// <param name="rightSide">The matrix to multiply.</param>
        /// <returns>The result of multiplication.</returns>
        /// <exception cref="ArgumentNullException">If <paramref name="leftSide"/> or <paramref name="rightSide"/> is <see langword="null" />.</exception>
        /// <exception cref="NotConformableException">If <strong>leftSide.Count != rightSide.Rows</strong>.</exception>
        public static DenseVector operator *(DenseVector leftSide, DenseMatrix rightSide)
        {
            if (leftSide == null)
            {
                throw new ArgumentNullException("leftSide");
            }

            if (rightSide == null)
            {
                throw new ArgumentNullException("rightSide");
            }

            if (leftSide.Count != rightSide.Rows)
            {
                throw new NotConformableException(Resources.ParametersNotConformable);
            }

            return (DenseVector)rightSide.LeftMultiply(leftSide);
        }
    }
}